Current noise of a superconducting single electron transistor coupled to a resonator 
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We analyze the current and zero-frequency current noise properties of a superconducting single 
electron transistor (SSET) coupled to a resonator, focusing on the regime where the SSET is operated 
in the vicinity of the Josephson quasiparticle resonance. We consider a range of coupling strengths 
and resonator frequencies to reflect the fact that in practice the system can be tuned to quite a 
high degree with the resonator formed either by a nanomechanical oscillator or a superconducting 
stripline fabricated in close proximity to the SSET. For very weak couplings the SSET acts on the 
resonator like an effective thermal bath. In this regime the current characteristics of the SSET are 
only weakly modified by the resonator. Using a mean field approach, we show that the current noise 
is nevertheless very sensitive to the correlations between the resonator and the SSET charge. For 
stronger couplings, the SSET can drive the resonator into limit cycle states where self-sustained 
oscillation occurs and we find that regions of well-defined bistability exist. Dynamical transitions 
into and out of the limit cycle state are marked by strong fluctuations in the resonator energy, but 
these fluctuations are suppressed within the limit cycle state. We find that the current noise of the 
SSET is strongly influenced by the fluctuations in the resonator energy and hence should provide a 
useful indicator of the resonator's dynamics. 

PACS numbers: 85.85.+j, 85.25.Hv, 73.50.Td 



I. INTRODUCTION 

The dynamics of a resonator coupled to a supercon- 
ducting single-electron transistor (SSET) is rather rich 
with a range of different behaviors expected to occur. 
Recent experiments in which the resonator was formed 
by a nanomechanical beam demonstrated ultra-sensitive 
displacement detection and cooling of the mechanical mo- 
tional 2 ". In another experiment^ the resonator was formed 
by a superconducting stripline and the SSET was used 
to drive it into a laser-like state 4 -^ of self-sustained os- 
cillation. 

A SSET consists of a superconducting island which is 
connected to two superconducting leads via tunnel junc- 
tions and capacitively coupled to a gate electrode, as 
shown schematically in Fig. [TJ The gate electrode forms 
part of the resonator, either as a metal wire deposited on 
top of a nanomechanical bearn^ or by forming the central 
electrode of a superconducting stripline^. Depending on 
the gate and bias voltages applied, the SSET can support 
a wide variety of current carrying processes^. Here we 
focus on a particular current resonance, the Josephson 
quasiparticle (JQP) cycl o 8 ' 9 ' 10 . At the JQP resonance 
current flows via a combination of the coherent tunneling 
of a Cooper pair at one of the junctions followed by the 
successive tunneling of two quasiparticles across the other 
junction. The center of the resonance occurs when the 
electrostatic energy of the two states linked by Cooper 
pair tunneling is the same. When the SSET is biased 
away from the center of the resonance the charges flow- 
ing through the SSET can either absorb energy from the 
resonato r 2 ! 11 ' 12 or emit energy into it&^ i 11 ! 13 ' 14 . Interest- 
ingly, the charge dynamics in the JQP cycle is closely 
related to that of another mesoscopic conductor namely 
a double quantum dot in the Coulomb blockade regime 
and so-called wide bias limi t 15 ' 16 . 



When the SSET is detuned from resonance in such 
a way as to emit energy into the resonator the latter 
is effectively pumped by the flow of charges. For suffi- 
ciently strong coupling the resonator can be driven into a 
state of self-sustained oscillation. Useful analogies can be 
made between the SSET-resonator system and quantum 
optical systems such as the laser— iliSil 7 .. In particular, 
there are a number of similarities between the predicted 
dynamics of a resonator driven by a SSET and a mi- 
cromaser. In a micromaser— a superconducting cavity 
interacts with a stream of two-level atoms which pass 
through it one at a time. The ordered flow of atoms 
which interact one at a time with the cavity leads to an 
interesting range of effects which are not seen in standard 
lasers such as the existence of a whole set of dynamical 
transitions and the regimes where the cavity is driven 
into non-classical state a 18 ' 19 . In the SSET-resonator sys- 
tem the resonator interacts with only one pair of charges 
moving through the SSET island at any one time, but 
the SSET current is not independent of the resonator 
dynamics. Nevertheless, features similar to the micro- 
maser such as the existence of a sequence of dynamical 
transitions and the possibility of driving the resonator 
into non-classical states are predicted to occur for the 
SSET-resonator system 4 ^ 14 ' 17 . 

In the micromaser system the state of the atoms emerg- 
ing from the cavity provides the information about the 
cavity dynamics 19 . Similarly, The SSET current pro- 
vides a natural source of information about the dynam- 
ics of the resonator—. However, as with many meso- 
scopic systems 2 ^, the fluctuations in the current contain 
much more information about the dynamics of the sys- 
tem than the average current alone. A number of recent 
studies 2 !^^^^^ have shown how the cur- 
rent noise of a nanoelectromechanical system can be used 
to infer quite a lot about the dynamics of the system. 
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FIG. 1: (a) Schematic diagram of the system consisting of a 
superconducting island formed by two tunnel junctions and a 
gate capacitor incorporating a resonator, (b) The JQP cycle 
involving coherent Cooper pair oscillations between the |0) 
and 1 2) charge states and incoherent quasi-particles tunneling 
to go from the |2) state to the |0) state via the |1) state. 



For example, it has been recognized that a bistability in 
the dynamics of the mechanical resonator can lead to an 
extremely large peak in the current nois o 25 ' 27 . In the 
case of the SSET-resonator system it has already been 
shown that the onset of self-sustained oscillations in the 
resonator can be associated with strong features in the 
SSET current noised. 

In this paper we present a systematic study of the cur- 
rent characteristics of a SSET when it is coupled to a 
resonator as a function of the SSET-resonator coupling 
strength, the choice of SSET bias point and the frequency 
of the resonator. In particular, we use a numerical ap- 
proach based on the master equation to explore the re- 
lation between the resonator's state (measured by the 
average resonator energy and associated variance) and 
the SSET current and zero-frequency current noise. For 
sufficiently strong coupling, strong correlations develop 
between fluctuations in the resonator energy and the cur- 
rent noise. This means that measurements of the current 
noise could provide a very useful probe of the resonator 
dynamics. 

In addition to the full numerical calculations, we also 
use a series of simpler approximate methods to gain 
more insight into the coupled dynamics of the system. 
When the SSET-resonator coupling is sufficiently weak, 
the SSET acts on the resonator like an effective thermal 
bath. In this regime we use a mean field approach that 
allows us to include information about the resonator dy- 
namics as well as the correlations between the SSET and 
the resonator progressively and hence discover how these 
affect the current noise without influencing the average 
current. In the strong coupling regime, where the res- 
onator undergoes oscillations driven by the current, we 
use an eigenfunction expansion of the Liouville operator 
in the master equatio n 24 i 31 ' 32 i 33 to understand the cur- 
rent noise. In the vicinity of a bistability the current 
noise is dominated by the slow switching of the system 
between the two effective states of the system which is 
manifested by one eigenvalue for the Liouvillian which is 
much smaller (in magnitude) than all the others 2 ^. In- 
terestingly, we find elsewhere that the noise can also be 
approximated quite well using a single term in the eigen- 
function expansion even when a wide separation between 
the smallest few eigenvalues does not exist. 



The organization of this paper is as follows. In Sec. 
HJwe introduce the master equation we use to model the 
SSET resonator system. We also describe how the steady 
state properties of the resonator and the current noise can 
be calculated numerically. Then in Sec. IIIII we present 
calculations of the SSET current and zero-frequency cur- 
rent noise together with the associated resonator energy 
and energy fluctuations for a wide range of system pa- 
rameters. We then focus on the weak coupling regime 
in Sec. IIVI where we present details of how simple mod- 
els based around the mean field equations of the system 
can be used to understand the current and noise in this 
regime. Then in Sec. |V] we consider the regime where 
the coupling is strong enough to drive the resonator into 
limit cycle states. We use eigenfunction expansions of 
the relevant Liouville operator to explore the extent to 
which the presence of a very slow timescale in the res- 
onator motion affects the current noise. Finally in Sec. 
IVII we draw our conclusions. Appendixes [ATlCl contain 
further details on certain aspects of the calculations. 



II. MASTER EQUATION FORMALISM 

In the vicinity of the JQP resonance 8 ' 34 ' 35 the SSET is- 
land is confined by charging effects to one of three charge 
states as shown in Fig. [TJ These states correspond to 
the presence on the island of no excess charges, |0), one 
Cooper pair |2), or one quasiparticle, |1). The master 
equation describing the SSET and resonator at the JQP 
resonance is given byiili, 



Pit) 



h 

C P {t) 



[He, p(t)] + C qpP {t) + C d p{t) 



(1) 



The first term describes the coherent evolution of the 
density matrix under the Hamiltonian H co while the sec- 
ond and third terms describe the dissipative effects of 
quasiparticle tunneling and the resonator's environment 
respectively. The SSET operators are defined in terms of 
the three accessible charge states, 

j»d = |0)<0|, pi = |1><1|, p 2 = |2)<2|, 
0=10X21, ?1 = |1)(2|, tfy = \0)(l\. (2) 

The Hamiltonian, H co , written in terms of these opera- 
tors takes the form, 

H co = AEp 2 -^-(c+J) + ^ + ~mn 2 x 2 



+ mil x s x (pi + 2p 2 ) , 



(3) 



where AE is the detuning from the JQP resonance, Ej 
is the Josephson energy and the resonator has frequency 
f2, mass m, momentum operator p and position operator 
x. The final term represents the linear coupling of the 
resonator to the charge on the SSET island. The length 
scale x s is the shift in the resonator position due to the 
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addition of a single electronic charge to the island. The 
coupling strength is conveniently expressed in terms of 

the dimensionless parameter k — — , where Vd s is the 
drain source voltage and e the electron charge. Note that 
here we use the language and notation appropriate for 
a nanomechanical resonator, but the Hamiltonian takes 
essentially the same form for a superconducting stripline 
resonator 3 . 

Quasiparticle decay at the right hand junction is de- 
scribed by the superoperator C qp , 

£q P p(t) = r(<?i +q2^p(t)(ql+q\) - -z {pi + p 2 , p(t)} , 

(4) 

where T is the quasiparticle tunneling rate and {•, •} is the 
anticommutato r 14 ! 36 . For simplicity, we have neglected 
both the differences between the rates for the two quasi- 
particle decay processes and the (weak) dependence of 
the rates on the position of the resonator—. The final 
term in Eq. ([1]) represents the damping of the resonator 
by its external environment: 

Cdp{t) = - 7eiCt ^ n0 ( next + ^ [x, [x, p{t)]] 

-^[xA^Pit)}}, (5) 

where j ex t is the damping rate and n ext = (e ti[1 / knText — 
where T ext is the temperature of the resonator's 
surroundings. 

The whole master equation can be represented by the 
single superoperator C, known as the Liouvillian. The 
Liouvillian operates in Liouvillc space where a Hilbcrt 
space operator a becomes a vector \a)) and both pre- 
(left) and post- (right) multiplication of the operator a 
can be represented by an appropriate matrix multiplying 
|a) ) 24 i 31 i 32 i 33 i 37 . The inner product for two vectors in 
Liouville space is defined as ((a\b)) = Tr[a^b\. Using this 
notation Eq. |T]) takes the form, 

\p(t)))=C\ P (t))). (6) 

Since we are dealing with an open system, the Liouvillian 
is not Hermitian and hence has different right and left 
eigenvectors, 

£M) = A p |r p », (7) 
((l p \C = X p ((l p \. (8) 

We choose to label the set of eigenvalues such that |Ao| < 
Ai| < .... Neglecting the possibility of degeneracy^ we 
assume that the eigenvectors form a complete orthonor- 
mal set, ((l p \r q )) = Tr[^r g ] = S pq . The solution to 
Eq. ^ can therefore be expanded in terms of the eigen- 
vectors to give 

Ip(*))> = E»(°)» eApt M> 

p=0 

= M+£((Zp|p(0))}e^|rp», (9) 
P =i 



where p(Q) is the initial density matrix of the system. 
For a master equation with a well-defined steady state 
(such as the one we consider here) the lowest eigenvalue 
will be Aq = 0, a property which we used to obtain the 
second line above. The other eigenvalues must obey 3 -^ 
Re (A p >o) < and the steady state density operator is 
\Pss}} — \p(°°))} = l r o))- The normalization of |ro)) is 
determined by Tr[p(i)] = 1, which gives ((Iq\ = where 
I is the identity operator (in Hilbert space). While |ro}} 
corresponds to the steady state, the eigenvectors \r p )) for 
p > each represent a change to the steady state density 
matrix that decays exponentially with rate — Re(A p ). 

The problem of finding the steady state density matrix 
is reduced to finding the right hand eigenvector of C cor- 
responding to the eigenvalue Ao = 0. By truncating the 
oscillator basis, Eq. ^ can be solved numerically to find 
the first few eigenvalues and eigenvectors of C Details of 
the numerical method and the approximations made are 
contained in Appendix [A"l 

Our aim in this paper is to understand to what ex- 
tent information about the dynamical state of the res- 
onator becomes imprinted on the transport properties 
of the SSET. As well as calculating the current we also 
consider the zero frequency current noise, which is in- 
dependent of the junction at which it is measured. We 
choose to calculate the noise at the junction at which the 
Cooper pair tunneling takes place (the left hand junc- 
tion) as the current operator here is composed of system 
operators alone, which along with the Markovian nature 
of the master equation, allows the use of the quantum 
regression theorem 1 ^. (The results have also been calcu- 
lated for the right hand junction and are in agreement.) 

The symmetrized current noise at the left hand junc- 
tion is defined as, 

Si l i l M = l°°dT (({I L (t + T),I L (t)})-2 (I L (t)) 2 ) , 

( 10 ) 

where the current operator at the left hand junction, i^, 
can be calculated by considering the flow of charge across 
the left hand junctional. The operator is given by 

Il = 1 ^(J-c). (11) 

For the current noise a symmetrized current operator can 
be defined in Liouville space, 

T L \p{t))) = \{I L p{t) + p{t)I L ) . (12) 

Using the quantum regression theorem to evaluate the 
correlation function, the current noise can be written in 
terms of Liouville space operators as 

S IlIl {u) =2[dT(((lo\I L e CM lL\ro)) - ((Zol^M) 2 ) , 

J OO 

In the zero frequency limit this has the solution 2 ! 

S lLlL (0) = -4((l Q \I L 1ll L \r )), (14) 
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where 7Z — QC~ 1 Q is the psuedoinverse of the Liouvillian 
and the projector Q = 1 — |?~o)) ((^o I- With this projection, 
the inversion takes place only in the space where C is 
regular. The current noise is conveniently parametrized 
by the Fano factor, which is defined as 



Fi 



2e(I) ■ 



(15) 



where 2e (I) is the Poissonian or shot noise limit cor- 
responding to the current noise for a single tunnel 
junction^. 



III. FEATURES OF THE CURRENT NOISE 

In this section we give a survey of the current and 
noise characteristics of the SSET and the corresponding 
resonator dynamics calculated numerically over a range 
of different parameters. In later sections we will provide 
a more detailed analysis of the most interesting regimes. 

The SSET-resonator system is rather complex, even 
within the framework of the simple model we are using 
here. In particular, the behavior of the system depends 
on a rather large number of different parameters. Some 
of the system parameters such as the detuning AE and 
the SSET-resonator coupling k can be changed during 
a particular experiment, while many of the other system 
parameters can be tuned by appropriate design of the de- 
vice, including the frequency of the resonator which can 
be in the region of 10MHz for a mechanical device^ or of 
order 10GHz for a superconducting stripline resonator 3 . 
We have not attempted a systematic survey of all feasi- 
ble parameter regimes, but instead focus primarily on the 
effects of changing AE, n and the resonator frequency. 

We start by reviewing the characteristics of the SSET 
in the uncoupled limit k — > 0. The current, (J)'" -0 , and 
Fano factor, Ff=°, for a SSET tuned to the JQP reso- 
nance are given by—'—: 



{iy 



4AE 2 
2 



h 2 T 2 



8E}> 



Ej 



3Ej> 

2h 2 r 2 



(AAE 2 + h 2 T 2 



3Ejf 



(16) 
(17) 



The current has a peak at the center of the resonance 
AE = 0, which has a width determined by T and Ej. 
Far from resonance the current Fano factor has a value 
of 2. This is because the rate at which the Cooper 
pairs tunnel onto the island is much slower than the 
quasiparticle decay rate and hence when a Cooper pair 
reaches the island it swiftly breaks up into quasiparti- 
cles. The two quasiparticle tunneling processes occur in 
quick succession^ (compared to the rate of Cooper pair 
tunneling) and hence the charge is effectively transferred 
in units of 2e. However, close to the center of the res- 
onance in the regime where Ej < hT (which we study 
here) there is a strong interplay between the coherent 





Z 
1.B 
1.6 
1.4 
1.Z 

'< 1 
D.a 
o.e 

0.4 
0.2 

-0.3 -D.15 D.15 

FIG. 2: (color online). Average energy of the resonator as a 
function of the detuning from resonance and coupling strength 
for fi/r = 0.12, Ej/eVd s = l/16,-y ext /T = 0.0001, r = 1 and 
next = 2. The dashed lines indicate transitions between dy- 
namical states: for most of the range considered the resonator 
is in the fixed point state, but for large enough coupling a 
transition to the limit cycle state occurs close to the center of 
the resonance. The bistable region is the smallest and occurs 
for k > 0.0011 and AE/eV ds ~ -0.15. 



transfer of Cooper pairs and the quasiparticle tunneling 
which results in a suppression of the noise. This suppres- 
sion is strongest at the center of the resonance where the 
coherent motion of Cooper pairs is most important. 

Coupling a resonator to the SSET can significantly 
modify the behavior of the two individual systems. In 
particular, energy can be transferred between the SSET 
charges and the resonator. For low to moderate cou- 
pling the resonator reaches one of three types of steady- 
state^^: a state in which the resonator fluctuates about 
a fixed point, a limit-cycle state where the resonator un- 
dergoes self-sustained oscillations and a 'bistable' state, 
where the two coexist. At larger couplings other states 
can be found for this system such as multiple limit cy- 
cles^, but these do not occur for the parameters used 
here. The different resonator states are readily identi- 
fied from the steady state number state distribution of 
the resonator, P(n) — Tr[\n)(n\p ss ], where \n) is a Fock 
state of the resonator. The fixed-point state has a single 
peak in the P(n) distribution at n = 0, while the limit 
cycle has a peak at n > and we define the bistable state 
as having two peaks. 

The behavior of the resonator and the correspond- 
ing current characteristics of the SSET for a slow res- 
onator, fi/T < 1, are illustrated in Figs. [2][5] The 
average energy of the resonator is shown in Fig. [5] as 
a function of the detuning AE and the coupling k for 
fi/r = 0.12. We have chosen a junction resistance 
r = Rje 2 /h = 1 and the quasiparticle decay rate is taken 
to be T = Vds/eRj. The Josephson energy is assumed to 
have a value Ej/eVd s = 1/16 so that throughout we will 
be in the regime where Ej < hT and hence the quasi- 
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FIG. 3: (colour online). Average current ((/) /eT) through 
the SSET as a function of the detuning from resonance and 
coupling strength. The dashed lines indicate transitions in the 
resonator's state. The parameters are the same as in Fig. [2] 
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FIG. 5: (color online). Current Fano factor, Fj, as a function 
of the detuning from resonance and coupling strength. The 
dashed lines indicate transitions in the resonator's state. The 
parameters are the same as in Fig. [2] and the colors are on a 
log scale. 
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FIG. 4: (color online). Resonator Fano factor, F n as a func- 
tion of the detuning from resonance and coupling strength. 
The dashed lines indicate transitions in the resonator's state. 
The parameters are the same as in Fig. and the colors are 
on a log scale. 



particles should be the dominant source of dephasing for 
the SSET island charge. The damping 'jext/F = 1 X 10~ 4 
(which is somewhat higher than might be expected in 
experiment) is chosen to ensure numerical convergence 
and a small amount of thermal noise has been included 

The transitions between the three different dynamical 
states of the resonator are indicated by dashed lines in 
Fig. [21 For AE < energy is transferred to the resonator 
and for strong enough coupling the resonator is driven 
into the limit cycle state which grows in size continuously 
as AE becomes more negative. However, for k > 0.0011 
when AE is sufficiently negative (AE/eVd s — —0.15) the 
resonator enters the bistable regime and then undergoes a 
transition back to the fixed point state in which the limit 



cycle disappears abruptly 4 -. The corresponding behavior 
of the current is shown in Fig. [3] Although the current is 
clearly modified by the coupling to the resonator, it does 
not contain any clear signatures of the transitions in the 
resonator state. 

An important measure of the resonator state is the 
Fano factor of the resonator occupation number, defined 
as F n — (An 2 )/(n) where (An 2 ) = (n 2 ) — (n) (where 
here n is the number operator a* a). F n is plotted in 
Fig. [3J Unlike the average energy of the resonator, F n is 
strongly peaked around the transitions between the fixed 
point and limit cycle states, with the strongest feature 
occurring in the vicinity of the bistable region 4 . 

The current Fano factor Fj is plotted in Fig. [SJ The 
behavior is rather complex, especially for relatively weak 
couplings, but overall it is clear that the current noise is 
a much better indicator of the presence of transitions in 
the resonator state than (I). The behavior is simplest 
for larger k, well within the regime where the dynamical 
transitions occur and in this region we see well defined 
peaks in the noise at the transitions into and out of the 
fixed point state. The noise peak is particularly promi- 
nent in the case of the bistability. Although we have 
defined the bistable state on the basis of just the number 
of peaks in the P(n) distribution, rather than the coexis- 
tence of two well separated states it is certainly possible 
to find regimes where a true bistability in this sense ex- 
ists (i.e. where the P(n) distribution not only has two 
peaks, but also has very small values for at least some 
range of the n values between the peaks). For a well de- 
fined bistability the current noise can become extremely 
large, a phenomenon which has been shown to be due to 
the existence of a very slow timescale associated with the 
switching of the system between the metastable states of 
the syste m 27 i 39 i 40 . 
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We now turn to consider the opposite regime of a fast 
resonator, fi/T ^ 1. In this regime the coherent coupling 
between the resonator and the SSET is expected to give 
rise to well-defined features when the resonator frequency 
matches an eigenenergy of the SSET, 



kttl = ±\ AE 



E 2 j 



±AE, 



(18) 



with fc an integer (for the relatively strong quasiparti- 
cle decay rates considered here fi/T ^> 1 means that 
Ej <C fi.fi). Numerically we do indeed find that limit 
cycles occur at these resonance points, almost always 
via continuous transitions, with the higher order features 
(|fc| > 1) appearing at progressively larger couplings. 

Figure [5] shows (I) , Fj and F n around the k = — 1 
(one photon absorption) peak corresponding to AE = 
— yj (HQ) 2 — E 2 j with the dashed line marking the onset 
of the limit cycle state. This resonance was recently 
observed in experiments using a superconducting res- 
onator—. It is clear that the peak in the current correlates 
well with the presence of the limit cycle, but this peak 
is present at the resonance even when the limit cycle is 
not formed, albeit with a very small size. However, the 
behavior of the noise shows a clear signature of the onset 
of the limit cycle with a distinctive peak structure form- 
ing along the transition lines for both Fj and F n (Fig. [6}d 
and c). Within the limit cycle regime both Fi and F n 
show dips the size of which is an indication of how well 
defined the limit cycle is. 

In practice it is also possible to build devices where the 
resistance of the junction where quasiparticle tunneling 
occurs is much larger than the quantum of resistance 3 
(i.e. r 3> 1). Increasing the resistance of the junc- 
tion where quasiparticle tunneling occurs enhances the 
coherent interaction between the Cooper pairs and the 
resonator. Strictly speaking one would expect other ef- 
fects which give rise to dephasing of the SSET charge 
(beyond just quasiparticle tunneling) to become relevant 
in this regime. Nevertheless, for our simple model we 
find that when r 3> 1 and fi/T ^> 1, the current noise 
at AE ~ — ^/(fifi) 2 — E 2 can become sub-Poissonian 
(Fi < 1) and the resonator Fano factor can also drop 
below unity, indicating a non-classical resonator state. 
Interestingly, this regime is quite distinct from the case 
discussed in Ref . |j (for which r ~ 1 and fi ~ T) where F n 
can also drop below unity. We note, however, that the 
thresholds for the sub-Poissonian regimes for F n and Fj 
are not perfectly correlated: these two effects can occur 
at the same time or separately depending on the param- 
eters chosen 2 ^. 



IV. THERMAL RESONATOR 

The simplest regime for the SSET-resonator system is 
that of very weak coupling where the resonator remains 
in the fixed point state for all values of the detuning 
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FIG. 6: (color online). (I), and F n for fi/T = 10, ~( BX t/T = 
0.0003 and n ex t = 0; the other parameters are the same as in 
Fig. [3 For the region within the dashed lines the resonator is 
in a limit cycle state and elsewhere it is in a fixed point state. 



AE. In this regime the effect of the SSET on the res- 
onator dynamics is analogous to an additional thermal 
bat h 11 ' 12 . In this section we consider in detail how the 
SSET current and noise are modified by the resonator in 
this regime and explore the extent to which the behav- 
ior can be understood in terms of simple models for the 
coupled SSET-resonator dynamics. 

For sufficiently weak coupling the resonator's steady 
state is a thermal (i.e. Gaussian) state to a good approx- 
imatio n 11 ' 12 . In this state the resonator position is deter- 
mined by the average charge on the SSET island which in 
turn is proportional to the average (steady state) current 
flowing through the SSET— ^ (see also Appendix |B| , 



W 2eT W 



(19) 



The average occupation number of the resonator, n, is 
given by a weighted sum of the contributions arising from 
the resonator's thermal environment, n ex ti and the effec- 
tive thermal bath it feels due to the SSET, n sse t, 



^Jext^ext i ^/sset^sset 
Jext + "Isset 



(20) 



The weighting factors are the resonator damping rates 
due to the true thermal bath, ^ ex t , and the effective bath 
due to the SSET, j sse t- For a slow resonator (fi <C T), 
both n sse t and 7 SS et are given by relatively simple ana- 
lytic expression o 11 ' 12 ' 41 : 



^fsset 



16mxjn 4 EjAE 



h 2 T 2 + AAE 2 
IQAEKl 



AAE 2 + l3h 2 T 2 + lOEj 
(AAE 2 + h 2 T 2 + 3Ejf _ 



(21) 
(22) 



Both jsset and n sset are strongly dependent on the de- 
tuning AE. Furthermore, 7 SS et becomes negative for 
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FIG. 7: (color online). Change in current through the SSET 
as a function of AE for k = 0.0001, Q/F = 0.05, Ej/eVds = 
1/16, "fext/F = 0.0001 and n ext = 2. The curves are labeled 
as num for the numerical results and fl for the change in 
the current calculated using Eq. (|24f) . [Note that for these 
parameters n varies from a value of 2 far from resonance to 
a peak value of 2.28 at AE = —0.01. r y S aet /jemt has maxima 
and minima of ±0.029 at AE = ±0.044.1 



AE < 0. However, for weak enough coupling the res- 
onator is stabilized in a thermal state by the damping 
arising from the coupling to the external bath. 

For weak SSET-resonator coupling the changes in the 
transport properties of the SSET due to the resonator 
are relatively small so it makes sense to examine just 
the difference between the values for the coupled and un- 
coupled cases. The change in the SSET current due to 
the resonator (calculated numerically) is shown in Fig. [7] 
We consider a slow resonator Q/T <C 1 and very weak 
coupling so that although the SSET has quite a strong 
influence on the resonator state, the resonator neverthe- 
less remains in a thermal state which is well described by 
Eqs. ([20 |) -([22 |) . From Fig.[7|we see that near the center of 
the resonance the current is suppressed by the resonator, 
but on either side of this there is an enhancement. The 
current noise is modified in a similar way to the current, 
but in the opposite sense, as shown in Fig. [8l thus there 
is an increase in the noise near to the resonance with a 
decrease on either side. 

The simplest way of including the influence of the res- 
onator on the SSET is to include the effect of fluctuations 
in the position of the resonator on the current. Because 
the resonator acts as a gate for the SSET island, a shift of 
the position of the resonator leads to an effective change 
in the detuning energy AE (Eq. [3]). Hence, when the 
resonator position fluctuates so will the detuning energy. 
We can incorporate the effect of the mechanical motion 
into the expressions for the current and noise, (Eqs. 1 1 64 - 
I17p , by calculating them for a fixed position, making the 
replacement AE — > AE + 2mVl 2 x s x, and then averag- 
ing over the resonator state. For the current Eq. HI 



fl 

- mcan3 




AE/cY^ 

FIG. 8: (color online). Change in the zero frequency Fano 
factor of the SSET due to the resonator. The curves are 
labeled as num for the numerical results, fl is obtained from 
Eq. (|26|l . mean2 is calculated using the second order mean 
field equations and meanS using the third order mean field 
equations. The parameters are the same as in Fig. [7] 



becomes: 
I(x) 



2eE 2 jT 



4 (AE + 2mQ 2 x s x) 



h 2 T 2 + 3E 2 



(23) 



Assuming the shift term is small we can perform a Taylor 
expansion and then take the average over the resonator. 
Keeping terms up to order x 2 s and using Eq. (|19p . we 
obtain 



\&mVL 2 x s AE . . 
1 - (x) 



16(mfl 2 x s 



fl 



P 

■ <*»>(! 



16A£^ 
~1~ 



24mQ 2 x 2 AE 



16(mn 2 x s ) 2 



(ieT 
(Ax 2 ) ( , 



(I) 



K, = 



16AE 2 
~1~ 



(24) 



where fl = 4:AE 2 + h 2 T 2 + 3E 2 and the averages are taken 
over the (Gaussian) steady state probability distribution 
for the resonator. The value of (Ax 2 ) is calculated using 
Eq. (HO]). 

For the current noise we naively replace AE — > AE ± 
Si(x) 



2mn 2 x s x to obtain 



2c 



21 (x)- 



16e£4r (E 2 + 2h 2 T 2 ) 



(^4 (AE + 2mfl 2 x s x) 2 + h 2 T 2 + 3E 2 

(25) 

After expanding to second order in x s and taking the 
average over the resonator state we can then calculate 
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the corresponding Fano factor, 



FP = 2- 



fj-48mn 2 x s AE 



32AE- 



(3-\QmVL 2 x s AE{x) - I6(mn 2 x s ) 2 (Ax 2 )(l 



16AE 2 



(26) 



where 4> = 8E 2 (E 2 + 2h 2 T 2 ) . 

Looking first at the current (Fig. [7]), it is clear that 
Eq. ([24]) accurately describes the modification due to the 
presence of the resonator. Thus in this weak coupling 
regime where the resonator remains in a thermal state, 
the modification of the current is simply due to the shift 
in the resonator's position (which gives the asymmetric 
shape) and a smearing out of the JQP current peak due 
to fluctuations in the resonator position. In contrast, we 
can see from Fig.[5]that for the Fano factor Eq. (f2"6"|) docs 
not capture the behavior correctly. Although the qual- 
itative shape is the same with a central peak with dips 
either side, the curves do not match and the asymmetry 
of the numerical curve is in the opposite direction to that 
predicted by the simple model. 

The reason for the disagreement in the current noise is 
that the simple model of a fluctuating gate neglects both 
the correlations between the electrical and mechanical 
motion and the dynamics of the resonator. The current 
noise (in contrast to the average current) is sensitive to 
the correlations between the SSET charge and the res- 
onator motion and hence to describe it accurately we 
need to include them in some way. A straightforward 
and systematic way to include correlations and informa- 
tion about the resonator dynamics can be found using the 
mean field equations of the system, namely the equations 
of motion for the expectation values of the SSET and res- 
onator operators. The mean field equations are generated 
in turn by multiplying the master equation by an oper- 
ator (or product of operators) and taking the trace over 
the full systemic. The mean field equations for the SSET 
resonator system are given explicitly in Appendix [5] up 
to second order. 

The set of mean field equations for the SSET res- 
onator system never forms a closed set with equations of 
motion for the operators always including some higher 
order operators. However, progress can be made by 
making a semiclassical approximation, in which corre- 
lations between certain operators are neglected so that 
the system of equations then becomes closed. Here 
we retain at least some of the SSET-resonator correla- 
tions by only applying the semiclassical approximation 
to products of higher order than required. For the set 
of second order equations we approximate products of 
three system operators, thus we make the substitution 

(x 2 c) -> 2 (x) (xc) + ((x 2 ) - 2 (a;) 2 ) (c). Crucially the 

correlations between products of any two operators are 
retained. 



The resulting set of equations is closed, but is non- 
linear because of the terms generated by the semiclassical 
approximation. However, by again using what we know 
about the steady state of the resonator in this regime [i.e. 
Eqs.[T9land[2"0] to replace the expectation values involv- 
ing only the resonator operators by their steady state val- 
ues, we can recover a linear set of equations (full details 
are given in Appendix |B| . This set of equations is then 
solved to obtain the steady state moments of the system 
using the same approach as that employed to solve the 
master equation (namely solving for the null eigenvector 
of the matrix of coefficients). The current is then ob- 
tained directly from the moments of the SSET operators, 
(I) = eT ((pi) + {p2})- The current noise can be calcu- 
lated using a slightly more elaborate calculation which 
introduces an electron counting variable^ 2 -, full details of 
which are given in Appendix [Cl 

The results from the second order mean field equations 
for the current agree very well with the numerical calcu- 
lation (and Eq. |2"1)) as one would expect. For the current 
noise the second order mean field calculation is quali- 
tatively correct as can be seen in Fig. [8l capturing the 
asymmetry in the numerical results though quantitative 
agreement is still lacking. This is not surprising as the 
second order mean field calculation only partly includes 
the SSET-resonator correlations and does not describe 
the dynamics fully. However, the mean field approach 
is readily extended to third order (i.e. the semiclassical 
approximation is only applied to products of four oper- 
ators), thereby including higher order correlations and 
more information about the resonator dynamics. We 
find that the third order calculation leads to quantita- 
tively correct results as shown in Fig. [5] However, we 
also note that reducing the coupling reduces the impor- 
tance of the higher order correlations which the second 
order mean field calculation neglects. Figure [5] provides a 
clear illustration of this as it shows that the second order 
calculation becomes accurate for low enough k. 



V. TIME SCALES OF THE SYSTEM 

In this section we discuss the signatures of the res- 
onator dynamics in the current and current noise when 
the coupling is strong enough to drive the resonator into 
limit cycle states. In particular we investigate how simple 
approximations to the noise based on the eigenfunction 
expansion in Liouville space (Eq. [5]) can be used to give 
insights into the connections between the fluctuations in 
the resonator state and the current noise. 

Figure [10] shows the current calculated numerically 
as a function of AE for three very different values of 
the resonator frequency. For low resonator frequencies 
(O/r = 0.12), the current is slightly suppressed at the 
center of the JQP resonance and enhanced further away. 
This is qualitatively the same as was seen for weak 
coupling in Sec. IIV1 even though now the coupling is 
larger so the resonator is driven into a limit cycle state. 
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FIG. 9: (color online). Change in the zero frequency Fano 
factor of the SSET due to the resonator for k — 5 x 10~ 6 . 
All other parameters and labeling of curves are the same as 
Fig.® 



- 
-- 




FIG. 10: (color online). Current as a function of AE for 
different resonator frequencies Q/T — 0.12, 1, 10. In each case 
the values of k and y ex t have been chosen to ensure that the 
system reaches the limit-cycle state for at least some values 
of AE whilst still remaining at low enough energies to allow 
a numerical calculation. For Q/T — 0.12, k = 0.0015 and 
Jext/T = 0.0001; for fi/r = I, k = 0.005 andj ext /T = 0.0008; 
and for Q/T = 10, k = 0.003 and j ext /r = 0.0003. The other 
parameters are the same throughout: Ej/eVds ~ 1/16, r = 
l,n ex t — 0. 



For the fl = T case small peak features are seen for 
AE < at points corresponding to the k = —1,-2 
and —3 resonances in Eq. (p~8|) . In the high frequency 
case (f2/r = 10), the current is greatly enhanced at the 
k = — 1 resonance and relatively unchanged elsewhere. 

Figs. [TTI fTSl show the current noise calculated numer- 
ically for the same parameters. For fl/T = 0.12 and 
fi/r = 1 the two peaks in the current noise correspond 
in both cases to a continuous transition from a fixed point 
state to a limit cycle at AE ~ and the presence of a 



region of bistability near the second (larger) peak in Fj . 
In between these two peaks the system is in a limit cycle 
state. For the Sl/F = 10 case the two peaks in Fi both 
correspond to continuous transitions (from fixed point to 
limit cycle state) with the resonator in a limit cycle state 
between the peaks. 



A. Bistability 

As discussed in Sec. IIII1 the current noise contains 
more information about the dynamical state of the res- 
onator than the current alone. For example, the presence 
of a dip in the noise between two strong peaks gives a 
clear indication that the resonator is actually in a limit 
cycle stated. The current noise can also tell us about the 
types of fluctuations present in the system, and the time 
scales over which these fluctuations decay. This analysis 
is particularly clear in the case of a bistability. 

Current noise peaks for bistable regions in nanoelec- 
tromechanical systems, such as the charge shuttle, have 
been studied extensivel y 25 i 27 ' 39 i 40 . The current charac- 
teristics of a conductor coupled to a truly bistable sys- 
tem (i.e. one with only two accessible internal states) 
can be described by a model specified in terms of four 
parameters: the (different) currents associated with the 
two states 7i, I2, and the switching rates between them of 
F12, r 2 i. The current and current noise for this two-state 
model take the simple for m 25 ' 27 



Sw(0) 



^2lh + ri2^2 

F21 + Ti2 
4 (A/ 2 ) 



21 



(27) 
(28) 



where (A/ 2 ) = r 2 ir 12 (Ii -I 2 ) 2 /( r 2i +T 12 ) 2 is the vari- 
ance in the current. 

This two-state model can be applied to a more complex 
system in a bistable regime if the two metastable states 
are well separated so that the switching rate between 
the states is much slower than the other relevant time- 
scales^! From Eq. {28]) we can see how slow switching 
rates between the two states can lead to a large value for 
the current noise in this regime. However, we also note 
that when the two metastable states give rise to very 
different currents the large variance that results can also 
make an important contribution to the current noise. 

For certain parameters the noise at the bistable tran- 
sition in our system is very well described by this two 
state model (with the two metastable states being the 
fixed point and limit cycle). In practice this means that 
a single set of the four parameters ^21, Ti2 can be 

found which allow us to fit the current and current noise 
to Eqs. (|2T|) and (J25J) respectively. Furthermore, the same 
parameters can be used to calculate higher cumulants of 
the current which can also be compared with numeri- 
cal results^ 5 -. The required parameters can be extracted 
as follows. The relative probabilities of the two states 
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r2i/(r 2 i + ri 2 ), T 12 / (r 2i + r 12 ) are obtained by inspec- 
tion of the steady state probability distribution P{n). 
Setting those elements of the steady state density matrix 
which correspond to just one of the two states to zero 
and recalculating the current then allows the currents 
ii and I2 to be obtained. Finally the sum of the rates 
Ti2 + r 2 i can be determined by comparing the current 
noise (calculated numerically) with Eq. [25] 

The two-state model can only be applied when a true 
bistability exists in the sense described in Sec. lIIII fi.e. the 
P(n) distribution for the resonator steady state should 
have two peaks with a vanishingly small probability for 
some range of n values in between) which we find gen- 
erally occurs for > T. Using the methods described 
above, we found that the two state approximation can 
be used to describe the current and current noise for the 
bistable state seen around AE/eVd s — —0.5 in Fig. IT51 
(n/r = 1), but not for the one around AE/eV ds ^ -0.12 
in Fig. [TT] (O/r = 0.12), where there is significant over- 
lap between the limit cycle and fixed point states. For 
the former case we obtained further confirmation that the 
two-state model could be applied by checking that the nu- 
merically calculated third cumulant agreed with that ob- 
tained using the two state model (an approach discussed 
in detail in Ref. 25]) and also by checking that the small- 
est (non-zero) eigenvalue of the Liouvillian matched up 
well with the total rate Ti2 + r 2 i (as we discuss below). 

In an ideal experiment one would be able to monitor 
the current with sufficient time resolution to observe the 
slow switching between two distinct values of the current 
directly. However, measuring the current noise as well 
as the average current in a region where the theory pre- 
dicts a bistability would provide convincing evidence if 
agreement was obtained. One could also make use of fur- 
ther generic predictions of the two-state model^l, such as 
the presence of a Lorentzian peak in the finite frequency 
current noise (at zero frequency) with a width given by 

T\2 + r 2 i. 



B. Eigenvector expansions 

In general, we cannot describe our system in the vicin- 
ity of the dynamical transitions by a simple two-state 
model. As we have seen, even where the transition in- 
volves a region of coexistence between the limit cycle and 
fixed point states the states may not be well enough sep- 
arated for a two-state model to apply. Near the contin- 
uous transitions between the limit cycle and fixed point 
states there are clearly not just two states involved. How- 
ever, one element of the two state model which might 
be expected to apply more widely is the emergence of a 
single very slow timescale which dominates the current 
noise. In the case of the continuous transition such a slow 
timescale might result from the vanishing effective damp- 
ing ("fsset +7ext) of the system at the transition. In what 
follows we use the eigenvector expansion of the Liouvil- 
lian to investigate the extent to which the current noise 



at each of the dynamical transitions can be described by 
a single slow process. 

We begin by rewriting Eq. (|14p for the current noise in 
terms of the eigenvectors and eigenvalues of the Liouvil- 
lian, 

S lLlL (0) = -Aj2^({lo\lL\r p ))((l p \I L \r )). (29) 

p—1 y 

This should be compared with a similar expansion for the 
variance in the current also for the left hand junction: 

(All) = (It) (II) 2 

= J2(( l o\ZL\r p }){{l p \l L \r )}. (30) 

The variance is given by a sum over the same matrix el- 
ements as the current noise but this time unmodified by 
the eigenvalues, X p . Each of the eigenvectors of the Li- 
ouvillian \r p )) describe a change to (or fluctuation away 
from) the steady state that decays with a purely exponen- 
tial rate — Re(A p ) [see Eq. [9|. Thus, the matrix element 
({Io\Tl\t p )) {{l p \lL\ro)) can be thought of as the variance 
in the current due to a fluctuation of type p. We then see 
that the current noise consists of a sum over the variances 
due to each type of fluctuation, each divided by the rate 
at which that fluctuation decays. 

It is clear from Eq. ([29]) that if |Ai| < |A 2 | then we 
could expect the current noise to be dominated by the 
first term, which corresponds to the slowest timescale 
in the system. This is in indeed what happens when 
the system has a well-defined bistability. In this case 
an obvious connection can be made with an appropriate 
two state model (i.e. Eq. [28)) , with the relevant eigen- 
value corresponding to the sum of the rates — Ai = 
Ti2 + r 2 i and the numerator gives the current variance, 
((lo\lL\ri))((h\X L \ro)) = (AI 2 ). More generally, it is not 
just a slow timescale that is important. For a single 
term in the eigenvector expansion to accurately describe 
the noise, the matrix element divided by the eigenvalue 
((l Q \l L \r p )) ({l p \I L \r )) / \ p for p = 1 must be much larger 
than for all p > 2. 

In Figs. [TT| - [T3l we compare the full current Fano factor 
with approximations using just the first term in Eq. (|29[) . 
The peaks at the transitions are described quite well by 
just the first term in the eigenvector expansion. Away 
from the peaks, however, we find that the noise is not 
captured by the approximation based on the first eigen- 
vector. It is particularly clear in Fig. [12] that something 
is missing from this approximation. The features that 
are simply due to the SSET alone are not captured, such 
as the dip at AE — and the Fano factor of 2 far from 
resonance. We can understand this better by considering 
the meaning of the eigenvectors and eigenvalues of the 
Liouvillia n 31 ' 43 . 

In the limit k — > the resonator-SSET system be- 
comes uncoupled and the eigenvectors and eigenvalues 
of the system can be expressed in terms of those of the 
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FIG. 11: (color online). Fi for « = 0.0015, fi/r = 0.12, 
Ej/eV ds = 1/16, 7e^t/r = 0.0001, r = 1 and n Mt = 0. The 
curve labeled num shows the numerical value of the noise, app 
is the approximate value of the noise using the first term in 
Eq. (|29[1 , and app5+Fj =0 is the first five terms plus the noise 
for a SSET alone [Eq. [32]. 




FIG. 13: (color online). Ft for k = 0.003, Q/T = 10, 
Ej/eV d s = 1/16, Tezt/r = 0.0003, r = 1 and n e:ct = 0. The 
curve labeled num shows the numerical value of the noise, app 
is the approximate value of the noise using the first term in 
Eq. (|29[1 , and app5+Fj =0 is the first five terms plus the noise 
for a SSET alone [Eq. [32]. 
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FIG. 12: (color online). Fj for « = 0.005, fi/r = 1, 
Ej/eV ds = 1/16, 7 e:Et /r = 0.0008, r = 1 and n ext = 0. The 
curve labeled num shows the numerical value of the noise, app 
is the approximate value of the noise using the first term in 
Eq. (|29[1 , and app5+Fj~° is the first five terms plus the noise 
for a SSET alone [Eq.l32]. 



individual subsystems, namely the SSET and the res- 
onator. When the resonator is decoupled from the SSET 
it still remains coupled to the external bath and its small- 
est (non-zero) eigenvalues are integer multiples*^ of j ex t- 
Thus the smallest of these eigenvalues corresponds to 
the energy relaxation rate of the resonator, — Jext> an d 
hence we can infer that the corresponding eigenvector de- 
scribes fluctuations in the resonator's energy (something 
which we will justify further below). There are also a 
set of eigenvectors (and corresponding eigenvalues) that 
describe fluctuations in the SSET charge state. In the 
uncoupled regime the current noise of the SSET can be 



obtained using Eq. (|29p . with the sum running over just 
the SSET eigenvalues, though we already know the result 
will be given by Eq. |T7|1 . 

For the coupled SSET-resonator system we can still 
identify the eigenvalues and eigenvectors as correspond- 
ing to one or other of the subsystems by looking at their 
behavior for large detunings (i.e. large \AE\) where the 
systems are effectively decoupled. The first few eigen- 
values, which correspond to the resonator, are shown 
(for the slow resonator case il/T = 0.12) in Fig. [T"4l as 
a function of AE. These first few eigenvalues indeed 
converge towards — Jemtt —Zjext, ■ ■ ■ for large detunings. 
Thus at least for large detunings the first eigenstate, |ri)), 
should therefore represent fluctuations which change the 
resonator energy. This can be confirmed by comparing 
the resonator variance to an expansion in terms of the 
eigenvectors, 



(An 2 ) =£<<Z |A/> p ))(<yA/>o)), 



(31) 



p=i 



where N\p(t))) = np(t) = a^ap{t). The full calculation 
of the energy variance is compared with approximations 
based on the first term in the eigenvector expansion in 
Fig. [151 K is clear that only the first eigenvector is needed 
to describe the energy fluctuations for large detunings as 
we expect. However, the approximation based on the first 
eigenvector also describes the energy fluctuations rather 
well at the peaks where the transitions occur, but not 
in between where the resonator is in a limit cycle state. 
However, Fig. [T5] also shows that we can describe the 
energy fluctuations throughout by using more terms in 
the eigenvector expansion. 

We are now in a position to understand why the cal- 
culation of the current noise using just the first term of 
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FIG. 14: (color online). The four smallest (non-zero) eigen- 
values as a function of the detuning, AE. The parameters 
are the same as in Fig. 1111 The eigenvalues differ from each 
other by more less than one order of magnitude throughout 
and converge towards integer multiples of j ext for large \AE\. 
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FIG. 15: (color online). Energy fluctuations of the resonator, 
F n , as a function of AE. The three curves show the full nu- 
merical calculation, num, and approximations using just the 
first term, app and the first five terms, app5, of the eigenfunc- 
tion expansion [Eq. [31] , respectively. The parameters are the 
same as in Fig. 1111 



the eigenvector expansion works as well as it does and 
to see how and why this can easily be improved upon. 
Comparing Figs. [11] and [15] it is clear that the single- 
eigenvector approximation to the current noise matches 
the numerical results well around the two peaks marking 
the transitions (between the fixed point and limit cycle 
states) where the first term in the eigenvector expansion 
also describes the energy fluctuations in the resonator ac- 
curately. The fact that the first term in the eigenvector 
expansion does not capture the current noise far from 
resonance is not surprising as it only describes fluctua- 
tions in the resonator state and does not include the fluc- 
tuations of the SSET degrees of freedom. We can easily 



obtain better agreement for large detunings by extending 
our approximation to include the contribution of the un- 
coupled SSET, Fj =0 . Better agreement within the limit 
cycle regime can be attained by using sufficient eigenvec- 
tors in our approximation to ensure that the fluctuations 
in the resonator energy are described accurately. Thus 
we arrive at our final approximate expression for the cur- 
rent noise 



Fj ~ Ff~ 



P =i 



((l \I L \r p 



X P e (I) 



(32) 



where m should be large enough so that the correspond- 
ing number of terms can be used to calculate (An 2 ) ac- 
curately [via Eq. [31] . In this case we find m — 5 is suf- 
ficient, and the current noise calculated this way agrees 
very well at almost all points as shown in Figs. rTTHT3l 
The one area where good agreement is still lacking us- 
ing Eq. [32] is within the limit cycle region for fi/T = 10, 
shown in Fig. 1131 This is because we have simply used 
the uncoupled contribution to the current noise arising 
from the SSET eigenvectors. In fact these SSET terms 
are strongly modified due to the resonant absorption of 
energy by the resonator from the Cooper pairs at this 
point. 

From these approximations it is clear that in the vicin- 
ity of the resonator transitions the current noise is largely 
determined by the slow fluctuations in the energy of the 
resonator. This is because the current depends in the 
first instance on the resonator position and hence on the 
latter's energy (as this is slowly changing compared to its 
period). Thus the current fluctuations depend strongly 
on the fluctuations of the resonator energy, rather than 
those of higher moments of the resonator. Thus when 
(An 2 ) depends on more than one eigenvector, the cur- 
rent noise does too. 

It is important to note that even in the regions where 
including just the first term in the eigenvector expan- 
sion describes the current noise fairly well this is not 
simply because the associated eigenvalue is very much 
smaller than all the others. We can see from Fig. [14] 
that (for these parameters) an overwhelming difference 
between the slowest two eigenvalues never develops and 
from Fig. 1161 that the relative size of the corresponding 
matrix elements is important in causing the first term in 
the eigenvector expansion to dominate. 



VI. CONCLUSIONS 

We have investigated the current and current noise of 
a SSET coupled to a resonator and how they relate to 
the latter's dynamics. The steady state properties of the 
system and the zero-frequency current noise are readily 
calculated using a numerical approach based on a repre- 
sentation of the master equation as a matrix equation in 
Liouville space. Overall we found that the current noise 
varies widely depending on the precise choice of SSET 
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FIG. 16: (color online). Matrix element corresponding to the 
two smallest magnitude eigenvalues, Ai , A2 . The parameters 
are the same as in Fig. 1111 Large differences between the 
magnitudes of the two matrix elements are visible in the same 
regions as the peaks in the corresponding plot of the current 
noise, Fig. [TT] 

bias point, the resonator frequency and the strength of 
the SSET-resonator coupling. For sufficiently strong cou- 
plings, the SSET current noise is strongly influenced by 
the fluctuations in the resonator energy. In particular, 
the resonator energy displays strong fluctuations in the 
regions where transitions between dynamical states oc- 
cur and this behavior is reproduced in the current noise. 
This means that measuring the current noise could pro- 
vide clear signatures of dynamical transitions in the res- 
onator. 

In addition to the full numerical calculations, we used 
a range of approximate methods to provide further in- 
sights into the coupled dynamics of the system. For very 
weak SSET-resonator couplings the SSET acts on the res- 
onator like an effective thermal bath. We found that in 
this regime mean field equations for the system operators 
provided a convenient way of establishing the importance 
of the SSET-resonator correlations and the resonator dy- 
namics in determining the SSET current and noise. For 
stronger couplings we used eigenfunction expansions of 
the Liouvillian matrix to demonstrate the strong connec- 
tion between the energy fluctuations in the resonator and 
the current noise. In many cases the current noise is well 
approximated by just a few terms in the eigenfunction 
expansion, but we found that it cannot be approximated 
accurately without including all of the eigenvectors that 
are needed to describe the energy fluctuations in the res- 
onator state. 
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APPENDIX A: NUMERICAL METHOD 

This appendix describes in more detail the numerical 
method used to solve the master equation and the ap- 
proximations that are made. To find the steady state of 
the system numerically the basis of the resonator must 
be truncated. External damping sets a limit on the res- 
onator energy. We therefore use a Fock state basis for 
the oscillator truncated to N states, where N is chosen 
to be large enough that the probability for the resonator 
to have an energy larger than fiflN is negligible. In Li- 
ouville space the density matrix is a vector and £ is a 
matrix with dimensions 97V 2 x 97V 2 . 

To obtain the steady state density matrix we need the 
eigenvector corresponding to the zero eigenvalue, or null 
eigenvector, of the Liouvillian. In Sec. [V]we also require 
the first few nonzero eigenvalues with the corresponding 
right and left eigenvectors. Due to the large dimension, 
complex nature and unsymmetric structure of the Liou- 
villian the inverse iteration method is use d 44 i 45 . Starting 
from a random vector \v))o the iteration is 

\v)) l+1 = (C-eiy'lvh 




where we have expanded \v))i in terms of the eigenstates 
of the Liouvillian and X is the identity in Liouville space. 
Repeating the iteration, the value of \v))i will converge 
to the eigenvector of C closest to e. To find the null 
eigenvector e is set to 10 -16 so that the matrix to be 
inverted is not singular. To find the other eigenvalues 
and eigenvectors we must start from an initial guess for 
the eigenvalue to be found. This guess can be updated 
every few iterations until the solution is found. The effi- 
ciency of the algorithm depends on how close e is to the 
exact eigenvalue in comparison to the next closest eigen- 
value. If the eigenvalue is known to high accuracy then 
convergence is found in a single iteration 4 ^. 

In order to use the largest possible value of TV we make 
two approximations. These approximations are based on 
the knowledge that if certain elements of the density ma- 
trix are known to be zero in the steady state then they 
can be omitted from the calculation by removing the ap- 
propriate rows and columns of the Liouvillian. We note 
that these terms must further be very rapidly decaying 
for the current noise not to be affected by their omission. 

The first approximation we make is to neglect the den- 
sity matrix elements corresponding to the q\ , qi , q\ and 
q\ operators. This is valid since there is no coherence 
between the pi state and the p or p 2 states so these 
elements must decay to zero in the steady state. This 
approximation reduces the dimensions of the Liouvillian 
to 57V 2 x 57V 2 . 
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The second approximation is made in terms of the os- 
cillator basis. The coherence between resonator Fock 
states which are widely separated in energy is small. El- 
ements of the oscillator density matrix far from the di- 
agonal can therefore be neglected. To check that this is 
a valid approximation the largest value on the last in- 
cluded diagonal of the resonator density matrix is found. 
So long as this value is below 10 -8 the results are found to 
be indistinguishable from the exact solution. After mak- 
ing these approximations the problem can be solved for 
N « 200 using an inverse iteration method on a desktop 
computer. The limiting factor is the memory required to 
calculate (C — el) 1 \v))i in the iteration procedure. 



APPENDIX B: MEAN FIELD EQUATIONS 

This appendix describes in detail the mean field equa- 
tion approach used in Sec. IIVI to solve the SSET- 
resonator system in the weak coupling limit. The mean 
field equations up to first order in the system operators 
were calculated in Ref. [14J. Here we go further and work 
to second order in the first instance, 



Ej 



(x) 
(v) 

(po) 

(Pi> 

(P2) 



= (v) 



(Bl) 



= -n 2 (x) - x s n 2 [( Pl ) + 2 (pa)] - lext (v) (B2) 



= ^((c)-(c t ))+ r (Pi) 

= -r<pi) + r<pa) 



(c) = 



(ct) = 



2h 



2mn 2 



(xc) 



h 2, 

2mfl 2 x 
+ i z — 



(xct) 



(B3) 
(B4) 
(B5) 



(B6) 



(B7) 



& ((xc) - (act)) + T (arpi) + (v Po ) (B8) 



2h 



(xPo) 

(xpi) = -T (xpi) + T (xp 2 ) + (vpi) 

(XP2) 



(B9) 



.Ej 



(xc) = 



(xct) 



1 IK ( {xc) ~ ^ ^ " r {xp2) + {vp2) (B10) 

A_£/ r\ , . Ej , . . . . . 



2mn 2 



i— — ) (xc^) 



(x 2 c) + (vc) 
Ej 



(Bll) 



h 2 

,2mfl 2 x. 



i( x Po) - ( x Pi)) 



h 



(upo) = i^- ((«c) - (v^)) + T (upi) - n 2 (xp ) 
(vpi) 

(VP2) 



(vc) = 



2h 

- lext (vpo) (B13) 

- (T + lext) (Vpi) + T (vp 2 ) - fl 2 (xpi) 

- x s n 2 (p x ) (B14) 

- ft 2 (spa) - 2x s fl 2 (pa) (B15) 

■ A£; r "\ / v n2/ v 
2 - Test I (uc) - il (xc) 

E 



i2mCl x= 



(vxc) 



(B16) 



(vet) 



««Po) - («P2)) 
r (aWcT) , 



(B17) 



where here the averages imply a trace over the SSET 
and resonator weighted by the density operator. 



The first thing to note is that Eqs. (|B1[) and (|B2|) can 
be used to obtain a simple approximation for the average 
resonator position. In the steady state we fincU^, (x) = 
—x s [(pi) + 2 (p 2 )], but from Eq. (|B4|I we see that (p 2 ) = 
(pi). Using the fact that the average current is related to 
the average charge on the SSET island, we then readily 
find that (x) = -Zx s (I) /(2eT). 

In order to obtain a closed set of mean field equa- 
tions at second order we need to make a semiclassical 
approximation to eliminate the third order terms (e.g. 
(x 2 ct) , (vxc)). This is done by setting the third order 
cumulant** to zero. In this context we apply the semi- 
classical approximation to products of three operators a, 
b and c by assuming that 



(abc) = (a) (be) + (b)(ac) + (c)(ab) - 2 (a)(6)(c) 



provided a, 6 and c all commute. Where the operators 
involved do not commute the expectation value should be 
symmetrized appropriately in order for the commutation 
relations to be preserved. 

Applying the semiclassical approximation to the term 
(x 2 c), in Eq. (|B11[) . we make the replacement (x 2 c) — » 

2(x)(xc) + {^(x 2 ) ~2(x) 2 ^j(c) and similarly for (zV) 

( B12 ) in Eq. (|BT2)) . The result ing approximate equations are 
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given by: 

(xc) 



~ l h 2 ) {XC) 

E 



,2ra0 2 



[(x 2 )-2(x) 2 )(c) (B18) 



(set) 



AE + AmQ 2 x s (x) T\ , ,. 
1 ft 2 J (XC ) 



+ i- 



2h 

. 2m0 2 



(x 2 )-2<z) 2 ) <ct).(B19) 



These equations can be linearized by treating (x) and 
(s 2 ) as parameters. To a good approximation the aver- 
age resonator position is given by Eq. (|19p in the steady 
state (with the current evaluated in the k — > limit) and 
as discussed in Sec. IIV1 within the weak coupling regime 
it is also possible to use Eq. (|20p to obtain (x 2 ). 

The other third order terms we need to consider are 
(vxc) and (xvcn, which arise in Eqs. (|B16[) and (|B17|) 
respectively. Since x and v do not commute we must 
first rewrite the expectation values in the following way 
before expansion so that the commutation relations are 
obeyed. 

(vxc) = - ((vx + xv) c) - i — (c). (B20) 
2 2m 

Performing the expansion as before we can make the re- 
placement, 

— ((vx + xv) c) — > (x) (vc) + (v) (xc) + — (c) (xv + vx) 



2 (x) (v) (c) . 



(B21) 



These equations are readily linearized by treating (x) as a 
parameter and using the fact that (xv) + (vx) = (v) = 
when the resonator is in a thermal steady state. The 
same procedure can be followed for the (xvc'') term to 
give 



(vc) 



AE + 2mn 2 x s (x) T \ . . 

-i ~ h ^-le X t) (VC) 

Ej 

+ i-^ ((vpo) - (vp 2 )) - x s fl 2 (c) 



(vrf) 



- n 2 (xc) 

AE + 2m,n 2 x s (x) T 

h 2 



(B22) 

lext ) (VC*) 



E 

n 2 (ict) , 



(B23) 



which completes the set of linearized equations. In order 
to solve this set of equations, we write the moments as a 



vector p so that the equations of motion can be rewritten 
in the form p = Ap where A is the matrix of coefficients. 
The steady state for the moments is then obtained from 
the null eigenvector of A. 

This method is readily extended to obtain the third 
order mean field equations by instead setting the fourth 
order cumulant to zero and following the same procedure. 



APPENDIX C: CURRENT NOISE 
CALCULATION IN THE MEAN FIELD 



This appendix describes the calculation of the current 
noise within the mean field model using the counting 
variable approac h 23 i 24 i 42 . We carry this out in terms of 
the second order equations, but the method is easily ex- 
tended to higher orders. First we write the master equa- 
tion in a modified form, where the number of quasiparti- 
cles, m, to have tunneled across the right-hand junction 
(since t = 0) are included, 



p(m,t) = -- [Hca, p(m,t)] + C d p(m,t) 
n 



+ r( gi + g 2 )p(m-M)(g| + gf) 
-^{pi+p 2 ,p(m,t)}. (CI) 



Mean field equations can be calculated for this mas- 
ter equation in the same manner as before, where we 
now include a subscript to indicate the number of quasi- 
particles to have tunneled. The majority of the equations 
are the same but with a subscript m, with averages now 
defined by (•),„ — Tr sys [-p(m, t)], with the trace taken 
over the system operators but not m. The following equa- 
tions have a modified form, 

(Po) m = «c) m - (ct) J + r ( Pl ) m _ 1 (C2) 



(Pi), 



- r (Pl)m + r (P2> ro _i 



(C3) 



E 

( x Po) m = ^ {(xc) m - (xc^J + T (a;pi) ro _ 1 

+ (vPo) m (C4) 

( x Pl)m = - T ( x Pl) m + r ( x P2) m -l + ( v Pl) m ( C5 ) 

E 

( V P0) m = ^ (( vc ) m ~ ( Uct ) m ) + r (^l)m-l 

- fl 2 (xp ) m - -/ext (vp ) m (C6) 

(«Pl)m = -( T + lext)(vpi) m +T(vp 2 ) m _ 1 

- n 2 (x Pl ) m - x s Q 2 ( Pl ) m . (C7) 

Due to the normalization condition, the total probability 
that m electrons have passed since t = is given by 



P(m,t) = (p )_ + (pi}_ + (p 2 ) r 



(C8) 
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The current noise is obtained using the MacDonald for- 
mula^: 



S(u) = 2uje 2 / dtsm(ut) 



dt 



^2 m 2 P(m,t) 



m=0 
\ 2 -i 



- J2mP(m,t) 



\m=0 



dt sin(ujt) 



2T (m{pi) + m(p 2 )) 



e e 2 



. . 2we 2 r r . 
2e (J) ; — [m(pi , iuj) — m(pi , — iu>) 



+ m(p 2 , iui) — m(p2, — itjj)] , 
where we have defined, 

oo 

m ( a ) =J2 m ^rn > 



(C9) 



(CIO) 



The equation of motion for m is found by multiplying the 
m-resolved equations of motion by m and performing the 
sum. For terms involving (a) m _ 1 we can use: 



m ( a )m-i = m ( a ) + («) 



(C13) 



m=0 



and the equation of motion for m is 



m = Am + y, 



(C14) 



where y is a vector containing the relevant inhomoge- 
neous terms. Laplace transforming and rearranging: 



rh{z) = — (zl — A) 1 y, 



(C15) 



where I is the identity. The final solution is written in 
terms of the vectors, 



and its Laplace transform 



m(a, z) = / die z m(a). 
Jo 



(Cll) 



To solve for the current noise we make use of the ma- 
trix notation introduced for the mean field equations in- 
troduced at the end of Appendix IB1 We start by defining 
the vector m, 



m 



(C12) 



k+ = (iuI-A) y, 

k = (—iuil — A) 1 y. 
The current noise is then given by, 



S{u) = 2e(J) 



(C16) 
(C17) 

(C18) 



+2e 2 F (k+ipx) + k-( Pl ) + k+( P2 ) + k-(p 2 )) , 

where k + (p\) indicates the element of k corresponding 
to pi. 
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